import warnings
notebook_warnings = []
def warning_collector(message, category, filename, lineno, file=None, line=None):
notebook_warnings.append(f"{category.__name__}: {message} (File {filename}, line {lineno})")
warnings.showwarning = warning_collector
# -*- coding: utf-8 -*-
# Copyright 2021 - 2022 United Kingdom Research and Innovation
# Copyright 2021 - 2022 The University of Manchester
# Copyright 2021 - 2022 Technical University of Denmark
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Authored by: Jakob S. Jørgensen (DTU)
# Edoardo Pasca (UKRI-STFC)
# Laura Murgatroyd (UKRI-STFC)
# Gemma Fardell (UKRI-STFC)
# Hannah Robarts (UKRI-STFC)
Key points demonstrated in this notebook:¶
Use CIL data readers to read in data¶
Use CIL Processors to manipulate, reduce and preprocess projection data¶
Use CIL Plugins for
ASTRAorTIGREtoolbox for forward and back-projection¶Use FBP for filtered back-projection reconstruction¶
Use CIL display tools
show2Dandislicerto visualise data and reconstructions¶Use iterative algorithms such as
SIRTas alternative if bad data¶Modify image geometry to reduce reconstruction volume to save memory and time¶
Data-set used in this notebook:¶
If you are running the notebook locally, install the data (usb.zip) from: https://zenodo.org/record/4822516#.YvJW5vfTXu0¶
First import all modules we will need:
import numpy as np
import os
import matplotlib.pyplot as plt
from cil.io import ZEISSDataReader, TIFFWriter
from cil.processors import TransmissionAbsorptionConverter, CentreOfRotationCorrector, Slicer
from cil.framework import AcquisitionData
from cil.plugins.astra import FBP
from cil.utilities.display import show2D, show1D, show_geometry
from cil.utilities.jupyter import islicer, link_islicer
Load the 3D cone-beam projection data of the USB:
# Please set the filename yourself, if you are running the notebook locally:
filename = "/mnt/share/materials/SIRF/Fully3D/CIL/Usb/gruppe 4_2014-03-20_1404_12/tomo-A/gruppe 4_tomo-A.txrm"
data = ZEISSDataReader(file_name=filename).read()
The data is loaded in as a CIL AcquisitionData object:
type(data)
cil.framework.acquisition_data.AcquisitionData
We can call print for the data to get some basic information:
print(data)
Number of dimensions: 3 Shape: (801, 1024, 1024) Axis labels: (<AcquisitionDimension.ANGLE: 'angle'>, <AcquisitionDimension.VERTICAL: 'vertical'>, <AcquisitionDimension.HORIZONTAL: 'horizontal'>)
Note how labels refer to the different dimensions. We infer that this data set contains 801 projections each size 1024x1024 pixels.
In addition to the data itself, AcquisitionData contains geometric metadata in an AcquisitionGeometry object in the geometry field, which can be printed for more detailed information:
print(data.geometry)
3D Cone-beam tomography System configuration: Source position: [ 0. , -95.04837036, 0. ] Rotation axis position: [0., 0., 0.] Rotation axis direction: [0., 0., 1.] Detector position: [ 0. , 55.08220291, 0. ] Detector direction x: [1., 0., 0.] Detector direction y: [0., 0., 1.] Panel configuration: Number of pixels: [1024 1024] Pixel size: [0.06585435 0.06585435] Pixel origin: bottom-left Channel configuration: Number of channels: 1 Acquisition description: Number of positions: 801 Angles 0-9 in radians: [-3.1415493, -3.1337836, -3.125914 , -3.1180005, -3.1102147, -3.1023476, -3.0945017, -3.0866468, -3.0787883, -3.0709336] Angles 791-800 in radians: [3.070867 , 3.0787256, 3.086563 , 3.0945282, 3.1023111, 3.110163 , 3.1180265, 3.1258702, 3.1337245, 3.1416595] Full angular array can be accessed with acquisition_data.geometry.angles Distances in units: units distance
CIL can illustrate the scan setup visually from the AcquisitionData geometry:
show_geometry(data.geometry)
<cil.utilities.display.show_geometry at 0x7fefb9548e30>
We can use the dimension labels to extract and display 2D slices of data, such as a single projection:
show2D(data, slice_list=('angle',220))
<cil.utilities.display.show2D at 0x7fefb8f9fce0>
From the background value of 1.0 we infer that the data is transmission data (it is known to be already centered and flat field corrected) so we just need to convert to absorption/apply the negative logarithm, which can be done using a CIL processor, which will handle small/large outliers:
data = TransmissionAbsorptionConverter()(data)
We again take a look at a slice of the data, now a vertical one to see the central slice sinogram after negative logarithm:
show2D(data, slice_list=('vertical',512))
<cil.utilities.display.show2D at 0x7fef9ad1e1e0>
Crop data by 200 pixels on both sides to save memory and computation time¶
data = Slicer(roi={'horizontal':(200,-200)})(data)
show2D(data, slice_list=('vertical',512))
New geometry: 3D Cone-beam tomography System configuration: Source position: [ 0. , -95.04837036, 0. ] Rotation axis position: [0., 0., 0.] Rotation axis direction: [0., 0., 1.] Detector position: [ 0. , 55.08220291, 0. ] Detector direction x: [1., 0., 0.] Detector direction y: [0., 0., 1.] Panel configuration: Number of pixels: [ 624 1024] Pixel size: [0.06585435 0.06585435] Pixel origin: bottom-left Channel configuration: Number of channels: 1 Acquisition description: Number of positions: 801 Angles 0-9 in radians: [-3.1415493, -3.1337836, -3.125914 , -3.1180005, -3.1102147, -3.1023476, -3.0945017, -3.0866468, -3.0787883, -3.0709336] Angles 791-800 in radians: [3.070867 , 3.0787256, 3.086563 , 3.0945282, 3.1023111, 3.110163 , 3.1180265, 3.1258702, 3.1337245, 3.1416595] Full angular array can be accessed with acquisition_data.geometry.angles Distances in units: units distance Shape out: (801, 1024, 624)
New geometry shape: (801, 1024, 624)
<cil.utilities.display.show2D at 0x7fefb8ef94c0>
CIL supports different back-ends for which data order conventions may differ. Here we use the FBP algorithm from the ASTRA Toolbox, which requires us to permute the data array into the right order:
data.dimension_labels
(<AcquisitionDimension.ANGLE: 'angle'>, <AcquisitionDimension.VERTICAL: 'vertical'>, <AcquisitionDimension.HORIZONTAL: 'horizontal'>)
data.reorder(order='astra')
data.dimension_labels
(<AcquisitionDimension.VERTICAL: 'vertical'>, <AcquisitionDimension.ANGLE: 'angle'>, <AcquisitionDimension.HORIZONTAL: 'horizontal'>)
The data is now ready for reconstruction. To set up the FBP algorithm we must specify the size/geometry of the reconstruction volume. Here we use the default one:
ig = data.geometry.get_ImageGeometry()
We can then create the FBP algorithm (really FDK since 3D cone-beam) from ASTRA running on the GPU and reconstruct the data:
fbp = FBP(ig, data.geometry)
recon = fbp(data)
show2D(recon, slice_list=[('vertical',512), ('horizontal_x', 325)], fix_range=(-0.1,1))
<cil.utilities.display.show2D at 0x7fef9a8904a0>
Offset initial angle to align reconstruction¶
# Note: The ZEISSDataReader reads the angles as radians, so we need to set the angle unit here:
data.geometry.set_angles(data.geometry.angles, initial_angle=-11.5*np.pi/180, angle_unit='radian')
<cil.framework.acquisition_geometry.AcquisitionGeometry at 0x7fef9a967200>
fbp = FBP(ig, data.geometry)
recon = fbp(data)
show2D(recon, slice_list=[('vertical',512), ('horizontal_x', 325)], fix_range=(-0.1,1))
<cil.utilities.display.show2D at 0x7fef9aa02ea0>
Use interactive islicer to flick through slices¶
islicer(recon,direction='vertical',size=10, minmax=(0.1,1))
islicer(recon,direction='horizontal_x',size=10, minmax=(0.1,1))
Extract and reconstruct only central 2D slice¶
data2d = data.get_slice(vertical='centre')
data2d.dimension_labels
(<AcquisitionDimension.ANGLE: 'angle'>, <AcquisitionDimension.HORIZONTAL: 'horizontal'>)
show2D(data2d)
<cil.utilities.display.show2D at 0x7fef9a606300>
ig2d = data2d.geometry.get_ImageGeometry()
recon2d = FBP(ig2d,data2d.geometry)(data2d)
show2D(recon2d,fix_range=(0,1))
<cil.utilities.display.show2D at 0x7fef9a6e2f00>
Simulate limited angle scenario with few projections¶
idx = [*range(50,400,10)] + [*range(450,800,10)]
# A number of other projection index ranges tried
# idx = [*range(50,350,10)] + [*range(450,750,10)] #+ [*range(400,500)] + [*range(600,700)]
# idx = [*range(50,150,10)] + [*range(200,350,10)] + [*range(450,550,10)] + [*range(600,750,10)]
# idx = [*range(25,375,10)] + [*range(425,775,10)]
# idx = [*range(0,125,10)] + [*range(275,525,10)] + [*range(675,800,10)]
# idx = [*range(50,200,5)] + [*range(300,450,10)] + [*range(550,800,20)]
# idx = [*range(100,350,10)] + [*range(500,750,10)]
# idx = [*range(0,100,20)] + [*range(350,500,20)] + [*range(750,800,20)]
plt.figure(figsize=(20,10)).set_facecolor('xkcd:white')
plt.subplot(1,2,1)
plt.plot(np.cos(data2d.geometry.angles), np.sin(data2d.geometry.angles),'.')
plt.axis('equal')
plt.title('All angles',fontsize=20)
plt.subplot(1,2,2)
plt.plot(np.cos(data2d.geometry.angles[idx]+(90-11.5)*np.pi/180), np.sin(data2d.geometry.angles[idx]+(90-11.5)*np.pi/180),'.')
plt.axis('equal')
plt.title('Limited and few angles',fontsize=20)
Text(0.5, 1.0, 'Limited and few angles')
Manually extract numpy array with selected projections only¶
data_array = data2d.as_array()[idx,:]
data_array.shape
(70, 624)
data2d.as_array().shape
(801, 624)
Create updated geometry with selected angles only¶
ag_reduced = data2d.geometry.copy()
ag_reduced.set_angles(ag_reduced.angles[idx], initial_angle=-11.5*np.pi/180, angle_unit='radian')
<cil.framework.acquisition_geometry.AcquisitionGeometry at 0x7fef9a6297f0>
Combine to new AcquisitionData with selected data only¶
data2d_reduced = AcquisitionData(data_array, geometry=ag_reduced)
Reconstruct by FBP¶
recon2d_reduced = FBP(ig2d,ag_reduced)(data2d_reduced)
show2D(recon2d_reduced, fix_range=(-0.1,1))
<cil.utilities.display.show2D at 0x7fef9a307170>
Try iterative SIRT reconstruction¶
Now set up the discrete linear inverse problem Ax = b and solve weighted least-squares problem using the SIRT algorithm:
from cil.plugins.astra.operators import ProjectionOperator
from cil.optimisation.algorithms import SIRT
A = ProjectionOperator(ig2d, ag_reduced, device="gpu")
Specify initial guess and initialise algorithm¶
x0 = ig2d.allocate(0.0)
mysirt = SIRT(initial=x0, operator=A, data=data2d_reduced)
Run a low number of iterations and inspect intermediate result¶
mysirt.run(10, verbose=1)
show2D(mysirt.solution, fix_range=(-0.1, 1))
<cil.utilities.display.show2D at 0x7fef9a43f980>
Run more iterations and inspect¶
mysirt.run(90, verbose=1)
show2D(mysirt.solution, fix_range=(-0.1, 1))
<cil.utilities.display.show2D at 0x7fef9a7ec620>
Run even more iterations for final SIRT reconstruction¶
mysirt.run(900, verbose=1)
show2D(mysirt.solution, fix_range=(-0.1, 1))
<cil.utilities.display.show2D at 0x7fef9a516d20>
Add non-negativity constraint using input lower=0.0¶
mysirt_lower0 = SIRT(initial=x0, operator=A, data=data2d_reduced, lower=0.0, update_objective_interval=50)
mysirt_lower0.run(1000, verbose=1)
Exception in thread Thread-4:
Traceback (most recent call last):
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/threading.py", line 1075, in _bootstrap_inner
self.run()
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/tqdm/_monitor.py", line 84, in run
instance.refresh(nolock=True)
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/tqdm/std.py", line 1347, in refresh
self.display()
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/tqdm/notebook.py", line 157, in display
pbar.value = self.n
^^^^^^^^^^
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/traitlets/traitlets.py", line 716, in __set__
self.set(obj, value)
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/traitlets/traitlets.py", line 706, in set
obj._notify_trait(self.name, old_value, new_value)
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/traitlets/traitlets.py", line 1513, in _notify_trait
self.notify_change(
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/ipywidgets/widgets/widget.py", line 700, in notify_change
self.send_state(key=name)
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/ipywidgets/widgets/widget.py", line 586, in send_state
self._send(msg, buffers=buffers)
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/ipywidgets/widgets/widget.py", line 825, in _send
self.comm.send(data=msg, buffers=buffers)
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/comm/base_comm.py", line 144, in send
self.publish_msg(
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/ipykernel/comm/comm.py", line 42, in publish_msg
parent=self.kernel.get_parent(),
^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/bgb37495/miniconda3/envs/cil_test_demos/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 783, in get_parent
return self._shell_parent.get()
^^^^^^^^^^^^^^^^^^^^^^^^
LookupError: <ContextVar name='shell_parent' at 0x7ff1eb743f60>
show2D(mysirt_lower0.solution, fix_range=(-0.1, 1))
<cil.utilities.display.show2D at 0x7feed8193ef0>
Compare all reduced data reconstructions in tighter colour range¶
show2D([recon2d_reduced, mysirt.solution, mysirt_lower0.solution], title=["FBP","SIRT","SIRT nonneg"], num_cols=3, fix_range=(-0.3,0.5))
<cil.utilities.display.show2D at 0x7feed8296360>
Compare horizontal line profiles¶
linenumy = 258
show1D([recon2d_reduced,mysirt.solution,mysirt_lower0.solution],
slice_list=[('horizontal_y',linenumy)],
dataset_labels=['fbp','unconstrained','constrained'],
line_colours=['black','blue','orange'],
line_styles=['dotted','dashed','solid'],
size=(12,8))
<cil.utilities.display.show1D at 0x7feed895d730>
Go back to full data FBP reconstruction, adjust reconstruction geometry to save time and memory¶
show2D(recon2d,fix_range=(-0.1,1))
<cil.utilities.display.show2D at 0x7feed823a150>
print(ig2d)
Number of channels: 1 channel_spacing: 1.0 voxel_num : x624,y624 voxel_size : x0.041692699432373054,y0.041692699432373054 center : x0,y0
Reduce the number of voxels¶
ig2d.voxel_num_x = 200
ig2d.voxel_num_y = 500
print(ig2d)
Number of channels: 1 channel_spacing: 1.0 voxel_num : x200,y500 voxel_size : x0.041692699432373054,y0.041692699432373054 center : x0,y0
recon2d = FBP(ig2d, data2d.geometry)(data2d)
show2D(recon2d,fix_range=(0,1))
<cil.utilities.display.show2D at 0x7fef9a530b90>
Centre the reconstruction volume around the sample:¶
ig2d.center_x = 30*ig2d.voxel_size_x
ig2d.center_y = -40*ig2d.voxel_size_y
print(ig2d)
Number of channels: 1 channel_spacing: 1.0 voxel_num : x200,y500 voxel_size : x0.041692699432373054,y0.041692699432373054 center : x1.2507809829711916,y-1.6677079772949222
recon2d = FBP(ig2d,data2d.geometry)(data2d)
show2D(recon2d,fix_range=(0,1))
<cil.utilities.display.show2D at 0x7feed82a44d0>
Further reduce the reconstruction volume¶
ig2d.voxel_num_x = 100
ig2d.voxel_num_y = 400
print(ig2d)
Number of channels: 1 channel_spacing: 1.0 voxel_num : x100,y400 voxel_size : x0.041692699432373054,y0.041692699432373054 center : x1.2507809829711916,y-1.6677079772949222
recon2d = FBP(ig2d,data2d.geometry)(data2d)
show2D(recon2d,fix_range=(0,1))
<cil.utilities.display.show2D at 0x7fef9a6b75c0>
print(ig2d)
Number of channels: 1 channel_spacing: 1.0 voxel_num : x100,y400 voxel_size : x0.041692699432373054,y0.041692699432373054 center : x1.2507809829711916,y-1.6677079772949222
Increase voxel size by a factor¶
ig2d.voxel_size_x = 4*ig2d.voxel_size_x
ig2d.voxel_size_y = 4*ig2d.voxel_size_y
print(ig2d)
Number of channels: 1 channel_spacing: 1.0 voxel_num : x100,y400 voxel_size : x0.16677079772949222,y0.16677079772949222 center : x1.2507809829711916,y-1.6677079772949222
recon2d = FBP(ig2d,data2d.geometry)(data2d)
show2D(recon2d,fix_range=(0,1))
<cil.utilities.display.show2D at 0x7feed89fbe00>
Reduce number of voxels by same factor¶
ig2d.voxel_num_x = 25
ig2d.voxel_num_y = 100
recon2d = FBP(ig2d,data2d.geometry)(data2d)
show2D(recon2d,fix_range=(0,1))
<cil.utilities.display.show2D at 0x7feed897d910>